import sys
sys.path.append('/home/aistudio/external-libraries')

import meshio
import matplotlib.pyplot as plt
import numpy as np
from scipy.spatial import Delaunay
from meshpy.tet import MeshInfo, build, Options

def plot_points_3d(points):
    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(points[:, 0], points[:, 1], points[:, 2], s=0.1)
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('Z')
    # 显示图形
    plt.axis("scaled") 
    plt.show()
    
def remove_duplicate_points(points):
    unique_points = np.unique(points, axis=0)
    return unique_points

'''----------------------------- Read the .msh file -----------------------------------'''
mesh = meshio.read('original_room_symmetry.msh', file_format='ansys')

# # Get the node coordinates and cell information
points0 = mesh.points
print("从fluent网格中读取节点数：", points0.shape[0])

points = remove_duplicate_points(points0)
print("删除重复点后剩余节点数：", points.shape[0])

# points 输入散点；furthest_site 为True时，计算最远点；incremental 为True时，允许增量添加点；qhull_options 为qhull参数，具体可参考qhull
tri = Delaunay(points, furthest_site=False, incremental=False, qhull_options=None)

mesh_info = MeshInfo()
mesh_info.set_points(tri.points)
mesh_info.set_facets(tri.convex_hull)
mesh_info.set_elements(tri.simplices)
mesh = build(mesh_info, options=Options("p"))

nodes = []
for i, n in enumerate(mesh.points):
    nodes.append(n)
nodes = np.array(nodes, dtype=int)

cells = []
for j, c in enumerate(mesh.elements):
    cells.append(c)
cells = np.array(cells)

print("\nDelaunay四面体化完成")
print("\nNew mesh infos:")
print("Nodes:", nodes.shape[0])
print("Cells:", cells.shape[0])

# 保存vtk文件
vtkFileName = "room_tet_remesh.vtk"
mesh.write_vtk(vtkFileName)

# 使用meshio读取VTK文件
mesh_vtk = meshio.read(vtkFileName)

# 将网格对象保存为Gmsh支持的格式（例如.msh文件）
mshFileName = "room_tet_remesh.msh"
meshio.write(mshFileName, mesh_vtk, file_format='gmsh22', binary=False)
print("网格已成功保存为vtk和msh文件:", vtkFileName, mshFileName)
